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It is shown that the semimetallic state of the two-dimensional honeycomb lattice with a point- 
like Fermi surface is unstable towards a canted antiferromagnetic insulator upon application of an 
in-plane magnetic field. This instability is already present at the mean-field level; the magnetic field 
shifts the up- and the down-spin cones in opposite directions thereby generating a finite density of 
states at the Fermi surface and a perfect nesting between the up- and the down-spin Fermi sheets. 
This perfect nesting triggers a canted antiferromagnetic insulating state. Our conclusions, based on 
mean-field arguments, are confirmed by auxiliary field projective quantum Monte Carlo methods on 
lattices up to 12 x 12 unit cells. 

PACS numbers: 71.10.Fd, 71.10.Hf, 71.27.+a, 71.30.+h, 73.43.Nq, 87.15. ak 



I. INTRODUCTION 

Graphenc, or the physics of electrons on the honey- 
comb lattice, has recently received tremendous attention 
due to its semimetallic nature with low-energy quasipar- 
ticles behaving as massless Dirac spinors; a recent review 
can be found in Ref. [T]. A crucial point is the stabil- 
ity of this semimetallic phase to particle-hole pairing. In 
particular, research activities have been devoted to the 
investigation of magnetic field induced transitions as a 
function of magnetic fields [2J [3J H] and electronic corre- 
lations. [5JG3I71IS] The vanishing density of states at the 
Fermi energy protects the semimetallic state against weak 
correlations. Notably, a finite critical value of the repul- 
sive Hubbard interaction U/t is required to destabilize 
the semimetallic state in favor of an antiferromagnetic 
Mott insulator. (SJ H [7] 

In this paper, we will argue that the semimetallic state 
is unstable against the application of an in-plane mag- 
netic field. The mechanism behind this instability can 
be understood already at the mean-field level. [9j [lOj [IT] 
The magnetic field generates a finite Fermi surface den- 
sity of states. A Stoner instability for arbitrarily small 
values of the Coulomb repulsion arises from the perfect 
nesting of the spin split Fermi sheets. This triggers an- 
tiferromagnetic order with staggered magnetization per- 
pendicular to the applied magnetic field and the opening 
of a charge gap. That the application of an in-plane 
magnetic field in the continuum limit facilitates a spon- 
taneous symmetry breaking has already been pointed out 
in Ref. [3]. The purpose of this paper is to show that 
those mean-field arguments indeed capture the correct 
physics, since exact quantum Monte Carlo (QMC) simu- 
lations on the honeycomb lattice compare favorably with 
those mean-field results. 

The outline of the paper is as follows. In Sec. |TT] the 
model Hamiltonian is introduced and its mean-field so- 
lution is discussed. The projective QMC method for 
ground state properties and numerical results are pre- 
sented in Sec. |III| The last section, Sec. |IV| contains the 



summary and the conclusions, also with respect to the 
experimental relevance of our findings. 

II. MODEL HAMILTONIAN AND MEAN-FIELD 
TREATMENT 

Our starting point is the Hubbard model on the hon- 
eycomb lattice shown in Fig. [T| 

(a) ^j\^k (b) p— * 















y ^ 





.r \k 



FIG. 1: (a) Real and (b) reciprocal space lattice vectors of the 
honeycomb lattice, ai = ao (0, 1, 0), &2 — ao (0, 1/2, y/3/2) 
and bi = 27r/a (0, 1, -l/V$), b 2 = 2-K/a (0, 0, 2/VS), with 
ao being the lattice constant. The unit cell with the orbitals 
a and b is indicated by the dashed diamond shape. Filled 
(empty) circles denote sites on the same sublattice. 

H = H a + Hjj + H B , 

i,r,cr 

H " = U E E (™M,T - 1/2) (n-, U - 1/2) , 

l—a,b i 

l—a,b i,tr 

The electron operator a]^ (&| CT ) creates an electron on 
the orbital a (b) in the unit cell i and the associated 
electron-density operator is n\ a — a\ a a- ua (b\ Jbi,a)i f° r 
I = a (b). Owing to the bipartite nature of the lattice, 
hopping with matrix element t occurs only between the 
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a and the b orbitals of unit cells related by lattice vec- 
tor r in the three directions {0,ai — &2, —3.2}- The on- 
site electron-electron repulsion is denoted by U > and 
p a = ±1 for g =T;I- I n t nc present case of half filling 
the chemical potential vanishes. In the following, we set 
(g/2) /xb = 1. We have included only a Zeeman coupling 
to the magnetic field. Hence, for comparison with exper- 
iments, we can only consider setups with magnetic field 
orientations parallel to the lattice plane since only in this 
case we can neglect the orbital coupling. 

The Hamiltonian H gives rise to two bands, 



A„(k) = p n 



-ik-r 



(2) 



with p n — ±1 for n = 1,2, respectively. At half-band 
filling the Fermi surface consists of two points, K,K' in 
Fig. [2] with density of states 



(3) 



Here, N corresponds to the lattice size and lineariza- 
tion of the dispersion relation around K and K' yields: 
p(ui) <x co/t 2 for w<t. 
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FIG. 3: (Anti-)ferromagnetic susceptibilities X(ajfm(Q) f° r 
the magnetic fields B/t = 0.01, 1.0. 



magnetic susceptibility tensor are given by 

xp M (q) = x;[~(q) + x2"(q), 



(•5) 



and correspond, respectively, to the ferromagnetic and 
the antiferromagnetic alignments of spins within the unit 
cell. Nesting occurs at q — Q — (0, 0) and leads to a 
logarithmic divergence of the antiferromagnetic mode at 
finite values of the magnetic field (Fig. |3j. 
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FIG. 2: Visualization of the nesting of spin-up and spin-down 
Fermi surfaces. In case of (a) B = the spin bands collapse 
onto each other, whereas for (b) B > the bands are shifted 
by virtue of the magnetic field leading to the nesting relation 
A llT (k) = -A 2 , i (k). 

Prior to examining the mean-field Hamiltonian we 
demonstrate the instability of the non-interacting sys- 
tem when turning on the magnetic field. Consider the 
transverse magnetic susceptibility tensor: 



x£ _(q) Xab_(q) 
xta (q) xtb (q) 



_ (x$y$ x + -(q) 
xt (q) (q) 



(4) 



Here xjy(q) = f Q P dr(S+(q, r)5 z 7(-q, 0)) and the in- 
dices I, I' = a,b label the sublattices. /3 corresponds to 
the inverse temperature and the spin raising and low- 



ering operators read S^(q) 



-iq-i 



4i ai d 



and 

Sa (q) = J2i e_lq '"i j."i,T' respectively. Similar defi- 
nitions hold for the b sublattice. The eigenvectors of the 



and 



Xafm (Q) 



^ V /(A», x (k))-/(A m , T (k)) 



In the above, 



k n,m=l 



(k)-A„ a (k) 



du> 



p{oj - B) + p{uj + B) 



to 



A„, CT (k) = A„(k)+p CT S 



(6b) 



(7) 



are single-particle states of Ho + Hb and /(w) = jt^j 
is the Fermi function. In the low-temperature limit, one 
can approximate the integral of Eq. (|6b|) to obtain 



Xafm (Q) 



UB)hx{^), \B\>0 
[constant, B = , 



(8) 
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FIG. 4: (a) Parallel magnetization my and (b) staggered mag- 
netization m± vs U and £>, obtained by numerically solving 
the gap equations, (111. 



where W corresponds to the bandwidth. Clearly, the 
divergence of the transverse susceptibility in the an- 
tiferromagnetic channel stems from the nesting prop- 
erty, Ai ! |(k) = — A2,j.(k). At zero magnetic field, this 
instability is cut off by the vanishing density of sates 
p(uS) oc ui/t 2 . At B > the low-energy density of states 
is finite thereby revealing the nesting instability. 

Given the above instability, the mean-field Hamilto- 
nian is derived by assuming the magnetization m to be 
alternating on the sublattices: m; = (0, m±(— l) 1 , my) 
with the index I = 0, 1 referring to the orbitals in the 
unit cell. That is, the magnetization m has a constant 
component my parallel to the field axis and a staggered 
component m± in the plane perpendicular to the field. 



To decouple the interaction part of the Hamiltonian 
we rewrite it as Hjj = — |?7X)i=a 6 Si w ^ n ke- 
ing the spin-1/2 operator, and the fluctuation term of 
S? ( = [nij + (Sj ; — m;)] 2 is omitted. With this ansatz, 
the mean-field Hamiltonian mixes the up- and the down- 



spin sectors thereby yielding four quasiparticle bands: 

2 



H 



MF 



k n,m— 1 



+ -UN(ml + ml), 



with 



#n, m (k) =p m y (A n (k) 



Here, p m = ±1, B cS - ^ ^ — 3 . 

Aj_ = ^Um±. Minimizing the free energy with respect 
to mil and m± yields the gap equations 



B 



B e ff) 



A,, 



Ai. 
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EE 

k n=i j (x n (h) + B cS y + ai 



u 

6N 



Ur p^ + B^+p^-B^ 

6 J y/rf + Ai ^ ; 



m H 



I 2 

N EE 



(A„(k) + Beff) 



AN 



k n=i J (\ n (k) + B cS y + ai 



dui 



(p(u> + B cS ) - p(lo - B cS ))uj 

, 



(lib) 



Our mean-field results are plotted in Figs . g ^a) , [f^c) , 
and|6^e). At zero magnetic field, we observed as a func- 
tion of U/t the expected transition from the semimetal- 
lic state (m = 0) at U/t < U c /t ss 3.3 to the antifcr- 
romagnctic Slater insulator characterized by \m±\ > 0. 
The semimetallic state at B = is characterized by the 
spin degenerate dispersion relation, A n . CT (k), as shown in 
Fig.[6ja). Ramping up the magnetic field lifts this degen- 
eracy thereby producing nested Fermi sheets of opposite 
spin indices. Hence, and irrespective of the magnitude of 
U <U C , energy can be gained by ordering the spins in a 
canted antiferromagnet. The energy gain corresponds to 
the ga p whi ch in the weak coupling limit, and by virtue 
of Eq. (11a), takes the form 



A, S We -3/2Up(B cff )_ 



(12) 



The dispersion relation of this canted antiferromagnetic 
state is plotted in Fig.^c). To compare at best with the 
QMC simulations we consider the quantity 



A\k : uj) = Im 

7T 



Gl a {k^)+G\ h {k,u)\ (13) 



with a finite broadening. As apparent, the features with 
dominant weight follow the dispersion relation Ai i |(k) 
and A2,|(k) and a gap at the Fermi level is apparent. Due 
to the transverse staggered moment, mixing between the 
up and the down dispersion relations occurs thereby gen- 
erating shadow features following the dispersion relations 
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FIG. 5: Staggered magnetization \J S + ~ (Q) /N (a) below and (b) above the critical interaction strength, (c) Finite-size 
extrapolation of S + ~ (Q) /N . In the semimetallic phase, U/t — 2, the data are consitent with the onset of transverse staggered 
order at finite magnetization M. For comparison, we have plotted the U/t = 5 data in the absence of a magnetic field. This 
value of the Hubbard interaction places us in the antiferromagnetic Mott insulating state. 



of Ai.j,(k) and A2,j.(k). The intensity of the shadow fea- 
tures tracks m± . As apparent from Fig. [3] the growth of 
m± as a function of the magnetic field is countered by 
the polarization of the spins along the magnetic field. It 
is interesting to note that irrespective of U/t the maxi- 
mal value of mj_ and hence of the magnetic field induced 
gap is at B = 1 corresponding to an energy scale match- 
ing the position of the Van Hove singularity in the non- 
interacting density of states. At this point a maximal 
amount of energy can be gained by the opening of the 
gap. 

Particle-hole symmetry can be exploited to map the 
repulsive Hubbard model onto an attractive Hubbard 
model where U < tunes the transition from a semimetal 
to an s-wave superconductor. |12j The external magnetic 
field driving the semimetal to a canted antiferromagnet in 
the positive-?/ model translates to doping the negative-// 
model which triggers a transition to a uniform superfluid 
phase with s-wave pairing, accordingly. |13j In the case 
of B — the self-consistent equation for the staggered 
magnetization, Eq. (11a), is identical to the BCS gap 



equation for the attractive- U case at half filling. 



III. PROJECTOR QUANTUM MONTE-CARLO 
METHOD 



To confirm our mean-field results, we have carried out 
projector auxiliary field QMC calculations. This PQMC 
algorithm is based on the equation 



reader to Ref. [14 . In this canonical approach, we fix 
the magnetization 



M 



(15) 



rather than the magnetic field (just as a magnetic field 
would induce a magnetization) . The magnetization M in 
the QMC algorithm corresponds to the mean-field mag- 
netization mil. N„ corresponds to the total number of 
electrons in the spin sector a. Furthermore, due to the 
particle-hole symmetry, which locks in the signs of the 
fermionic determinants in both spin sectors one can avoid 
the so-called negative sign problem irrespective of the 
choice of the magnetization. In practice, for each finite 
system, we choose a value of the projection parameter 9 
large enough so as to guarantee convergence within sta- 
tistical uncertainty. 

To detect transverse staggered magnetic order under 
an applied magnetic field, we have computed the spin- 
spin correlation functions 

S + ~(q)= ^E e " lq(i ~ j) (^ + (i)^a(j)-^ + (i)^(j)> • 



i.J 



(16) 

If long-range staggered magnetic order perpendicular to 
the applied field direction is present, then 



S+ (Q) QMC 

lim = 77i i 

N^oo N - 1 - 



(17) 



(goMlgc 
<*o|*o) 



lim 



T \e 



-6H 



Ae 



-OH 



vfr,. 



(\Me" 



-28H\\f, 



(14) 



The trial wave function \iPt) has to be non-orthogonal 
to the ground state wave function, (iprl^Po) an d the 
ground state is assumed to be non-degenerate. For the 
details of the formulation of this approach, we refer the 



acquires a finite value. We have computed this quantity 
on 6 x 6, 9 x 9, and 12 x 12 lattices, and our results 
are plotted in Fig. [5] both for U < U c and U > U c . 
At U/t = 2 < U c /t and zero magnetization, M = 0, 
our results are consistent with = whereas, at 

M = 2/6, m2 MC takes a finite value. Although we can- 
not reproduce the essential singularity of the mean-field 
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FIG. 6: Single-particle spectral function ^4J(k, lo) at U/t = 2, based on the mean-field (left) and the QMC (right) calculations, 
respectively. The magnetization M takes the values of 0, 2/6, and 4/6 (from top to bottom). For the intermediate value of 
the magnetization, one can clearly see the opening of a quasiparticle gap and the avoided level crossing as m± acquires a finite 
value. In general, the magnitude of the gap tracks m± [Fig. [HJa)]. The QMC spectral functions were obtained via analytical 
continuation of the Green's functions with the stochastic maximum entropy method. For the QMC calculations, the lattice 
size was set to 12 x 12 unit cells. The legends at the top indicate the false color values. 
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U/t = 2,M = 2/6 



U/t = 2,M = 4/6 
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FIG. 7: The plot directly compares the QMC (red) spec- 
tral intensity profiles (k, w) against the mean-field solution 
(blue, dotted). Same parameters as in Fig. [6] 



calculation at U < U c , the overall form of the trans- 
verse staggered magnetization compares favorably with 
the mean-field results [see Figs. [5^ a) and [5jb)] both at 
U < U c and U > U c . 

Within the PQMC, the zero-temperature single- 
particle Green's function along the imaginary time axis 
can be computed efficiently with methods introduced in 
Ref. [15]. From this quantity, we can obtain the spectral 
function of Eq. (13) with the use of a stochastic formu- 



lation of the maximum entropy method. |16l 117] The so 
obtained results for ^(k,w) are plotted in Figs. [(3jb), 
[6jd), and [6^f). As apparent the features in the QMC 
calculation which are associated with substantial spec- 
tral weight are well reproduced by the mean-field calcu- 
lation. The particle-hole transformation, at 



and b, 



-b\ -a, leads to the relation 



,4 T (k,w) = A L (k,-uj). 



a-, -„ 



(18) 



At finite magnetic fields or equivalently at finite magneti- 
zations, the staggered transverse order leads to a gapless 
Goldstone mode of which the quasiparticle can spin-flip 
scatter. As a consequence, and as already observed in the 
mean-field calculation, the features of the down spectral 
function should be visible in (k, ui) . Upon inspection of 
Fig.[7]one will observe that for each dominant low-energy 
peak at w(k) in Al(k, us) a shadow feature at — w(k) is 
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FIG. 8: Charge gap Ax as a function of the applied magnetic 
field for t — 2.5eV, obtained by numerically solving the gap 
equations, (111. 



IV. 



SUMMARY AND CONCLUSIONS 



In conclusion, we have carried out mean-field calcula- 
tions and projective quantum Monte Carlo simulations 
for the Hubbard model on the honeycomb lattice in a 
magnetic field oriented parallel to the lattice plane. With 
this setup, only the Zceman spin coupling is present. Our 
results show the inherent instability of the semimetallic 
state to a canted antiferromagnet upon application of the 
magnetic field. Ramping up the magnetic field generates 
nested up and down Fermi surfaces with finite density 
of states. As a consequence, the onset of canted anti- 
ferromagnetic order opens a charge gap and provides an 
energy gain irrespective of the magnitude of the Hubbard 
repulsion. 

Experimentally, such a transition could be observed 
by magneto resistance measurements. The transition to 
the canted antiferromagnet breaks a U(l) symmetry and 
hence occurs at finite temperatures Tkt in terms of a 
Kosterlitz-Thouless transition. Below Tkt the power- 
law decay of the transverse spin-spin correlation func- 
tion should suffice to produce a visible pseudo-gap in the 
charge sector and hence an increase in the resistivity as 
a function of decreasing temperature. The primary is- 
sue to observe the transition is the magnitude of the re- 
quired magnetic field so as to obtain a visible gap. With 
t w 2.5 - 3.0eV and U « 10 - 16eV, [H] one can read- 
ily see that very large magnetic fields will be required 
to obtain charge gaps in the meV region. In particular, 
in Fig. [8] we plot the charge gap in meV as a function 
of B in tesla for values of the Coulomb repulsion close 
to U c . The g factor has been set to the (approximate) 
free-electron value, g = 2. As apparent, depending on U, 
values on the order of B oc 10 2 — 10 3 T are required to 
obtain an acceptable gap. Clearly those numbers imply 
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that the only feasible manner to observe this effect would 
be to grow graphene directly on a magnetic substrate. 

With those numbers in mind, we can now consider 
magneto resistance experiments carried out in layered 
highly oriented pyrolitic graphite. 3 For magnetic fields 
perpendicular to the plane and at low temperatures, an 
insulating state as characterized by ^ is observed at 
B oc 0.1T. For fields parallel to the plane, intensities of 
roughly B oc 10T are required to observe the insulating 
behavior in the magneto resistance. Given this data, it is 
clear that the dominant effect of the magnetic field stems 
from the orbital coupling rather than from the Zeeman 
spin coupling. [l|j] The authors of Ref. [3] account for the 
parallel field data by mentioning deviations from perfect 
alignment between the graphene plane and the magnetic 



field. Given our estimate of the required magnetic field to 
achieve a charge gap for the parallel field configuration, 
we can only confirm this point of view. 
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